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Abstract 



There have been suggestions that heat caused by cerebral metabolic 
activity may constrain mammalian brain evolution, architecture, and 
function. This article investigates physical limits on brain wiring and 
corresponding changes in brain temperature that are imposed by ther- 
modynamics of heat balance determined mainly by Na + /K + -ATPase, 
cerebral blood flow, and heat conduction. It is found that even mod- 
erate firing rates cause significant intracellular Na + build-up, and the 
ATP consumption rate associated with pumping out these ions grows 
nonlinearly with frequency. Surprisingly, the power dissipated by the 
Na + /K + pump depends biphasically on frequency, which can lead to 
the biphasic dependence of brain temperature on frequency as well. 
Both the total power of sodium pumps and brain temperature diverge 
for very small fiber diameters, indicating that too thin fibers are not 
beneficial for thermal balance. For very small brains blood flow is 
not a sufficient cooling mechanism deep in the brain. The theoreti- 
cal lower bound on fiber diameter above which brain temperature is 
in the operational regime is strongly frequency dependent but finite 
due to synaptic depression. For normal neurophysiological conditions 
this bound is at least an order of magnitude smaller than average 
values of empirical fiber diameters, suggesting that neuroanatomy of 
the mammalian brains operates in the thermodynamically safe regime. 
Analytical formulas presented can be used to estimate average firing 
rates in mammals, and relate their changes to changes in brain tem- 
perature, which can have important practical applications. In general, 
activity in larger brains is found to be slower than in smaller brains. 

Keywords: Metabolism; Heat balance; Brain size; Fiber diameter; Temperature; 
Wiring; Limits. 
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1 Introduction 



Brain like any computational device (Landauer 1961; Bennett, 1982) or biological or- 
gan (Rolfe and Brown, 1997) dissipates energy (Siesjo, 1978; Ames, 2000; Clarke and 
Sokoloff, 1994). Recent studies on cerebral metabolism indicate that brain is energy 
expensive (Aiello and Wheeler, 1995; Attwell and Laughlin, 2001; Lennie, 2003), as it 
consumes relatively more energy than most other organs, and its total metabolic rate 
has larger scaling exponent than whole body metabolic rate (Karbowski, 2007). Poten- 
tial thermal imbalance between heat produced and removed could lead to overheating 
and aberrant functioning, because neural properties are strongly temperature depen- 
dent (Koch, 1998). Therefore, the problem of heat transfer could be a serious factor 
shaping brain evolution (Falk, 1990), organization, and function (Baker, 1982; Raichle, 
2003; Kiyatkin, 2007). Consequently, what is the relationship between brain tempera- 
ture, its size, cerebral power generated and neural activity? Can we estimate changes 
in firing rates based on changes in temperature? Is there any thermodynamic limit on 
brain size? If so, does 5 kg, which is the mass of the largest mammalian brain (Haug, 
1987; Purves, 1988), approach that limit? Are neural sizes and fibers also constrained 
by thermodynamics, and if so to what extent? How fast can neural computation be to 
maintain thermal balance and physiological range of intracellular ionic concentrations? 

This article answers these questions by finding the power generated by sodium 
metabolic pumps and relating this power to the thermal and neuroanatomical prop- 
erties of brain tissue. From this, we determine theoretical thermal bounds on fiber 
dimensions, and compare them with empirical data. These bounds enable us to esti- 
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mate the upper theoretical limits on the density of brain wiring. 



2 Methods 



2.1 Voltage and Na + -K + dynamics. 



Na + and K + are two major ions affecting neural membrane dynamics (Kandel et al, 
1991; Koch, 1998) and its energetics (Astrup et al, 1981; Erecinska and Silver, 1989; 
Rolfe and Brown, 1997; Ames, 2000). The energy consuming Na + /K + pump affects 
membrane electrical properties, because it pumps out 3 Na + and pumps in 2 K + per 1 
ATP molecule consumed (Kandel et al, 1991). The equations governing membrane and 
Na + /K + -ATP activities of a neuron are given by: 



CS— = -g Na S(V - V Na ) - g K S(V 



V K ) - 9lS(V - V L ) -I p -l 




9NaS(V - V Na ) - 3I P - al. 



g K S(V - V K ) + 2I P - PL 



(1) 



where V is the membrane potential, C is the membrane capacitance per unit area, 
S is neuron's membrane surface area (primary axons and dendrites), V^ a and Vk are 

4 



Na + and K + reversal potentials, gNa and gx are volt age- dependent Na + and K + con- 
ductances per unit membrane area, ql is the leak conductance and Vl is the reversal 
potential corresponding to the leak current. The parameter U n denotes neuron's vol- 
ume, F is the Faraday constant. The current I p is the Na + /K + pump current given 
by I p = AS[Na] fc /([Na] fc + 9 k ) where 9 = 20.0 mM, A is the maximal pump current 
per membrane surface area, and k is the Hill constant. The symbol [...] denotes intra- 
cellular ionic concentration, I s is the synaptic current, and a, (3 are voltage dependent 
proportionality parameters determined below. 

Voltage dependent conductances are represented by g^ a = m 3 hg Na and gx = n A g K , 
where g Na and g K are the maximal Na + and K + conductances. The gating variables h 
and n obey the standard kinetic equation dx/dt = a x (l — x) — f3 x x, (for x — h, n), and 
the fast variable m is set to its equilibrium value a m / (cx m +/3 m ). Voltage dependences of 
the parameters a x and (5 X were chosen, with a slight modification, as in the Traub-Miles 
model (1991), which describes a pyramidal neuron in the hippocampus. Specifically, 
for sodium channels: a m = 0.32(V + 54)/ (1 - exp(-(V + 54)/4)), (3 m = 0.28(V + 
27)/ (exp((y + 27)/5) - 1), a h = 0.128 exp(-(V + 50)/18), (3 h = 4.0/(1 + exp(-(V + 
27)/5), and for potassium channels a n = 0.0AA(V + 52)/ (1 - exp(-(V + 52)/5)), f3 n = 
0.5exp(— (V + 57)/40), where V is expressed in mV. 

It is assumed that the neuron represented by Eq. (1) is embedded in the network of 
neurons firing with an average firing rate /. The synaptic current of a single synapse 
is qg s e~ l l Ts (y — V s ), where q is the release probability, g s is the maximal synaptic con- 
ductance (it is assumed that the rising phase of synaptic conductance is much faster 
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that its decaying phase that is characterized by the time constant t s ), and V s is the 
reversal potential for synapses. Since the majority of synapses in the mammalian brain 
are excitatory and glutamate (Braitenberg and Schuz, 1998), we have V s = (Koch, 
1998). The average synaptic conductance is / Jq^ dt qg s e~ t l Ts ! , which in the physiolog- 
ically valid limit fr s <C 1 yields qg s fT s . Neurophysiological data indicate (Markram et 
al, 1997) that the release probability q is frequency dependent. Here, we follow Dayan 
and Abbott (2001) and assume that synaptic depression modifies the release probability 
according to q(f) = qo/(l + 7^/), where go is the release probability at Hz, is the 
depression time constant, and 7 is the parameter controlling the degree of depression 
(for 7 = lack of depression, for 7 = 1 maximal depression). Thus, the total aver- 
age synaptic current is I s = fq Mg s T s (V — V s )/(1 + 7^/), where M is the number 
of synapses (or presynaptic neurons) per neuron. Fluctuations around the average I s 
cause the neuron to fire an action potential with an average firing rate /. Glutamate 
synapses with non-NMDA receptors are used more frequently for regular transmission 
(NMDA receptors are blocked by Mg 2+ ; Kandel et al, 1991) and therefore only these are 
considered here. Synaptic current of non-NMDA type is composed almost exclusively 
of Na + and K + ions, which is the reason for the presence of al s and f3I s terms in the 
dynamics of [Na] and [K] in Eq. (1), where a + (5 — 1. 

The parameter a can be determined in a few steps. First, we can write the synaptic 
current I s as I s = (jNaiY — Vno) + 9k(V — Vk), where g K and g^a are K + and Na + 
conductances (voltage independent) at glutamatergic synapses. Since for V — V s — 
the synaptic current I s = 0, we have c/k = —gNaVNal^K- Second, from the condition 
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al s = g Na (V - V Na ), we obtain a = g Na (V - V Na )/(g Na (V - V Na ) + g K (V - V K )), 
which leads to a = Vk(V — Vno) / {V{Vk — Vno,)}- Note that a is voltage dependent. 
Consequently, the synaptic contribution to changes in intracellular Na + concentration 
is given by al s = q(f)fMg s r s (V — V Na )V K / (V K — V Na ). The synaptic contribution to 
changes in intracellular K + is (3I S = -q(f)f Mg s r s (V - V K )V Na /(V K - V Na ). 

Equation for voltage dynamics in Eq. (1) assumes that the membrane is equipo- 
tential, i.e., axons and dendrites have on average equal potentials. This assumption 
neglects spatio-temporal effects associated with action potential propagation. The va- 
lidity of this approximation is discussed in the Discussion section and in Appendix 
A. 

Adding appropriately sides in Eq. (1) we obtain the "charge conservation" equation: 
UnF^ ([Na] + [K]) = CS^- + g L (V - V L ). (2) 



This equation reflects the fact that non-zero membrane potential is caused by concen- 
tration gradients of charged ions across membrane, and temporal changes in potential 
are directly related to the temporal changes in ionic concentrations. This temporal 
relationship during an initial phase of an action potential is discussed in the Results 
section. Here, we focus instead on the long-term relation between intracellular [Na] 
and [K]. In this regime and for firing rates < 100 Hz (duration of action potential and 
voltage recovery ~ 10 — 15 msec), the average of V is close to its resting potential V Q , 
which in turn is very close to Vl- Thus, on the long-time scale, the right hand side of 
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Eq. (2) is close to 0, which implies that [Na] + [K] const. That is, changes in [Na] 
directly determine changes in [K], and practically there is no need to solve the equation 
for the potassium dynamics. Values of neurophysiological parameters used in this paper 
are presented in Table 1. 

2.2 Neuroanatomical relationships. 

In studying metabolic and thermodynamic properties of brain tissue the following neu- 
roanatomical relations were used: (i) volume density of synapses p s = NM/U g is brain 
size independent, where N is the number of neurons and U g is the cortical gray matter 
volume (Braitenberg and Schuz, 1998; DeFelipe et al, 2002); (ii) the fraction 1 — <f> of 
the gray matter volume taken by fibers (axons and dendrites) is approximately constant 
and about 2/3 (Braitenberg and Schuz, 1998); (iii) white matter volume U w scales with 
gray matter volume U g as U w = O.lQQUg- 23 (Zhang and Sejnowski, 2000), where U g and 
U w are expressed in cm 3 . 

The relationship between neuron's surface area S and gray matter volume U g can 
be determined as follows. The gray matter volume is U g = NU n + (pU g , where neuron's 
volume U n = (l a d 2 a + Iddfyn/A with l a , l d denoting average (unmyelinated) axon and 
dendrite length in the gray matter and d a , dd denoting their corresponding diameters 
(neuron's soma is neglected as it is small). The parameter is the fraction of volume 
taken by non-fibers (synapses, blood vessels, etc). The neuron's surface area is S = 
Tt{l a d a + Iddd)- We define "an effective" fiber diameter d as 

d = (l a d 2 a + l d d 2 d )/(l a d a + l d d d ), (3) 
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which is proportional to the neuron's volume to surface ratio, i.e., d = 4U n /S. In fact, if 
volumes of axons and dendrites are equal, which seems to be neuroanatomically valid in 
the gray matter (Braitenberg and Schiiz, 1998; Chklovskii et al, 2002), then d represents 
a harmonic mean of axonal d a and dendritic dj diameters, that is d — 2d a dd/ (d a + dj). 
This formula can be derived by noting that axon length l a = ld(dd/d a ) 2 (equal volumes 
of axons and dendrites), and substituting l a to Eq. (3). Note that for constant d d) 
the effective diameter d decreases with decreasing d a . In the limit d a /dd <C 1, we have 
d w 2d a . 

For neuroanatomical values of d a = 0.3 /im and dd = 0.9 /im in mouse (Braitenberg 
and Schiiz, 1998), we obtain the effective fiber diameter d = 0.45 /im. Since the average 
values of cortical axon and dendrite diameters are roughly invariant with respect to brain 
size (i.e., they should be more or less the same for mouse and human; Braitenberg and 
Schiiz, 1998), we can expect that the effective fiber diameter d should not change with 
brain size either. Therefore, we will use the above value of d — 0.45 /im for illustrative 
purposes throughout the paper. Apart from that, occasionally we will also show results 
for the smallest experimentally known axons of diameter 0.1 /im corresponding to d — 
0.18 /xm (for d d = 0.9 /im), and for the smallest molecular ly feasible axons (neurites) of 
diameter 0.06 /im (axons without Natrium channels that lack action potentials; Faisal 
et al, 2005), corresponding to d — 0.11 /xm (also for dd = 0.9 /im). 

With the definition of d we can relate the total surface area of neurons, NS, to the 
gray matter volume U g as 
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NS = 4(1 - <f))U g /d, 



(4) 



and relate the surface density of synapses along axons and dendrites, Mj S, to synaptic 
density p s as 

M _ p s d 

-J ~ 4(1^)- (5j 

These relations are used below to rewrite Eq. (1) in a more convenient form for numer- 
ical and analytical calculations. 

2.3 Explicit voltage and ionic dynamics. 

With the modifications and interdependencies between parameters described in Sub- 
sections 2.1 and 2.2 we can rewrite Eq. (1) as: 



° d ~L = ~ 9Naiy ~ VNa) ~ 9k{y ~ Vk) ~ 9l{y ~ Vl) 

A[Na] fe q(f)fp s dg s r s V 

+ synaptic jluct. 



[Na] fe + 9 k 4(1 - 0) 
dFd[M ( q(f)fp s dg s r s V K \ 3A[Na] fe 

" - \9Na + -Tf. 77777 {V - VNa) ~ 



4 dt Y ±{l-4>){V K -V Na )) y iyaj [Ne] k + e k 

dFd[K] ( q(f)fp s dg s T s V Na \ 2A[Na] fc 



where the release probability modulated by synaptic depression is q(f) = qo/(l + 7t^/) 
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(Dayan and Abbott, 2001). These equations are solved both numerically (using the 
4th-order Runge-Kutta method) and approximately analytically. Fast synaptic fluctu- 
ations denoted symbolically in the top equation of Eq. (6) cause the neuron to fire 
stochastically with the average frequency / (i.e. on average, every 1/f seconds the 
voltage is set to —20 mV). 

2.4 Electric power generated by Na + /K + -ATP pumps. 

Since the metabolic rate in white matter is 3-4 times lower than in gray matter (Siesjo, 
1978; Karbowski, 2007), the white matter contribution to the total brain metabolic 
power is neglected. The average power Patp dissipated by iV neurons in the gray 
matter is given by 



where the integral represents the electrical work performed by Na + /K + -ATP pumps 
during the long time At, much larger than the average interspike interval 1/f. This 
work goes for removing 3 Na + ions and importing 2 K + ions against their electrochemical 
gradients, which cost 1 ATP molecule (Kandel et al, 1991). Opposite signs in front of 
3I p and 2I p indicate the fact that Na + and K + ions move in opposite directions through 
the membrane. 

Integral in Eq. (7) can be estimated by noting that in the long-time limit the pump 
current I p assumes its average value (with some fluctuations around, see the Results) 




(7) 
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given by I PAV = A5'[Na]^/([Na]^ + 9 k ), where [Na] at , is the average long-term sodium 
concentration (determined in the Results section). In this limit, the average value of 
voltage V(t) is approximately equal to its resting value V a . Thus, the electric power 
Patp is 



Patp ~ ([NaL + 6)fc) ^ ~ 2V « ~ V <>) > (») 



where = (RT/F) ln([Na]ex/[Na] „) and V K = (RT/F) \n([K] ex /[K] av ), with the 
extracellular Na + and K + concentrations (Hille, 2001): [Na] e:r = 145 mM, and [K] ex = 4 
mM. 

2.5 Mechanisms of brain cooling. 

The major contribution to the cerebral metabolic rate or heat constitute Na + /K + - 
ATPase (Astrup et al, 1981; Erecinska and Silver, 1989; Rolfe and Brown, 1997; Ames, 
2000). Therefore, in this paper, the heat coming from other reactions such as glycolysis 
(conversion of glucose to ATP) is assumed to be less important and hence it is neglected 
in the estimates. In deep brain regions heat generated by hydrolysis of ATP molecules 
and coupled to it activity of Na + /K + pumps is transfered by conduction and circulating 
cerebral blood flow, whereas at the surface (scalp) the heat is removed mostly by con- 
vection/conduction and radiation (Fig. 1). Heat conduction (convection) is associated 
with the existence of a temperature gradient between the brain or scalp and the external 
environment. Heat radiation is a quantum mechanical property of every physical object 
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having temperature above absolute zero. Transfer of heat to the circulating cerebral 
blood in deep brain regions is possible because incoming blood (artery) has a slightly 
lower temperature than the brain tissue in this region (Hayward and Baker, 1968; Nybo 
et al, 2002; Kiyatkin, 2007). In general, the incoming blood temperature corresponds 
to the core body temperature, which is assumed to be a constant because it changes on 
time scales that are much slower than changes in brain temperature (Kiyatkin, 2007). 

The equation governing heat balance and spatial distribution of brain temperature 
T(r) has the form (Nelson and Nunneley, 1998; van Leeuwen et al, 2000; Sukstanskii 
and Yablonskiy, 2006): 



where k is the arithmetic mean of the thermal conductance of brain tissue (including 
cerebrospinal fluid, skull, and scalp at the edge), p b \ and c b \ denote blood's density 
and specific heat, CBF is the cerebral blood flow rate expressed in sec -1 , T b \ is the 
incoming blood (arterial) temperature equivalent to body core temperature, and Gatp 
is the heat (per time unit) generated in the gray matter due to Na + /K + -ATPase. In 
general, Gatp > Patp because Patp is the useful work (per time unit) performed due 
to hydrolysis of ATP and release of Gatp of free energy (per time unit). In Sec. 3.5 we 
estimate the relative magnitude of Gatp and Patp, i-e., we calculate the efficiency of 
the sodium pump. The parameter U br denotes brain volume, i.e., U br = U g + U w (brain 
geometry is modelled as half of a ball; Fig. 1). The parameters p br and c br denote 




2 T 



p bl c bl CBF(T - T bl ) + GATp/Ubr, 



(9) 



dr 2 
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brain tissue density and specific heat and they are approximately equal to pu and cu, 
respectively. Only steady-state regime of Eq. (9), i.e. dT/dt = 0, is considered below. 
The boundary condition imposed on Eq. (9) has the following form: 



dT 



(asBiT^-T^+viT^-T^), 



(10) 



dr 



where R is the brain's radius given by R = (3U br /2n) 1 / 3 , a SB is the Stefan-Boltzmann 
constant, rj is the heat convection/conduction coefficient between the scalp and the 
outside environment, and T sc and T Q denote the scalp and environment temperatures. 
The left hand side of Eq. (10) corresponds to the heat rate removed from the scalp per 
surface area. The right hand side of Eq. (10) is the sum of scalp radiation (term ~ T S 4 C ) 
and scalp convection/conduction (term ~ T sc ). The presence of the term proportional 
to T% is a consequence of the fact that scalp not only radiates energy to the environment, 
but it also receives some radiation from it. 

The heat removed from the brain through the conduction is given by 
Q c = — k J d 3 r d 2 T/dr 2 = —2irR 2 KdT/dr\ r=R . The heat removed by the circulating 
cerebral blood is given by Qu = puCuCBF J d 3 r (T(r) — T u ). Therefore, by integrating 
Eq. (9) for the whole brain volume, at the steady-state, we obtain: 



GaTP — Qc + Qbl, 



(11) 
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where the conduction term Q c is a sum of scalp convection/conduction and scalp radi- 
ation, i.e., Q c = Q cv + Q r , with Q cv = 27rR 2 r)(T sc - T Q ) and Q r = 2irR 2 a S B(Tg C - T 4 ). 
Values of all thermodynamic parameters are presented in Table 1. 



3 Results 

3.1 Sodium influx during a single action potential. 

When synaptic fluctuations cause the neuron to fire an action potential, the Na + con- 
centration first rise and then slowly decays to its equilibrium value, due to the workings 
of the Na + /K + pump. Sodium influx during a single action potential can be estimated 
using Eqs. (2) and (6). It is composed of the two contributions (Fig. 2A): Na + influx 
during a rising phase of sodium conductance g^ a and voltage V, and Na + influx during 
a decline phase of g^ a and V. The total influx is given by (see Appendix A): 

A[Na] c » 4{C + d 6C \ v Na - V ), (12) 

where SC denotes "correction" to the effective capacitance coming from the prolonged 
Na + channels activation and is given by 5C ~ 0.064^ A r a r o (VAr a — 0.6Vk)/(Vno, — V a ). 
Value of r is taken from numerical simulations of Eq. (6), and is approximately r Q pa 0.4 
msec. This number enables us to determine the parameter SC in Eq. (12), which we find 
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to be 8C 2.4 /iF/cm 2 for typical resting values of voltages: Vn u — 68 mV, Vk = —100 
mV (for [Na]= 12 mM, [K]= 155 mM; Hille, 2001) and V Q = -67 mV. It should be 
remembered, however, that 5C is not a constant, but it varies slightly depending on the 
level of intracellular sodium concentration [Na] via V Na and V K (see Fig. 2B). 

The immediate conclusion from Eq. (12) is that the amplitude of the sodium influx 
increases with decreasing the effective fiber diameter d, and this agrees with a direct 
numerical integration of Eq. (6), as is shown in Fig. 2C. Based on the above values, we 
find from Eq. (12) that for d = 0.45 /im (harmonic mean of average axon and dendrites 
diameters in mouse; Sec. 2.2) the sodium influx A[Na] D ^ 0.42 mM, which corresponds 
to ~ 2.5 • 10 5 Na + ions per /im 3 . For a mouse neuron with equal volume of axons and 
dendrites, and with l a — 4 cm and d a = 0.3 /im, this gives the total influx of 14.2 • 10 s 
Na + ions. Direct numerical integration of Eq. (6) for d = 0.45 /im yields a similar 
sodium influx, A[Na] D ^ 0.47 mM, which corresponds to the total influx of 15.9 • 10 8 
Na + ions for a mouse neuron. 

The influx of 14.2 — 15.9 • 10 8 Na + ions during an isolated action potential obtained 
above is comparable to the estimate of Attwell and Laughlin (2001), who used a dif- 
ferent, phenomenological, approach and obtained a slightly lower value of 11.5 • 10 s 
Na + ions. The difference can be attributed to the differences in the assumption regard- 
ing spatial properties of the membrane potential and the amplitude of depolarization 
during an action potential. In this paper, it is assumed that the membrane is equipoten- 
tial, i.e., axons and dendrites are depolarized by the same amount (135 mV), whereas 
Attwell and Laughlin (2001) assume that dendrites are 50 % less polarized than ax- 
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ons (dendrites 50 mV, axons 100 mV). Below, it is shown that the formula for Na + 
charge influx used by these authors is equivalent to the formula (12), if we assume that 
the axon depolarization AV a during an action potential is the same as the dendrite 
depolarization AV d . The total charge influx Q AL in Attwell and Laughlin (2001) is 
Qal = 4ffC(l a d a A V a + IdddAVd), if we neglect a small soma contribution. The prefac- 
tor of 4 was chosen by these authors to a large extent arbitrary, and it comes from the 
effect of simultaneous activation of Na + and K + channels, which is analogous to the 
presence of the 5C contribution in Eq. (12). Now, assuming that AVd = AV a = AV, 
and noting that Qal = FU n A[N&]AL, where A[Na]^ L is the sodium influx in the Attwell 
and Laughlin (2001) formulation, we obtain AfNaj^L = 16CAV/(Fd), where d is the 
effective fiber diameter defined in Eq. (3). The formula for AfNaj^L is very similar to 
the formula (12), except for the numerical factor in front, which is 13.6 in Eq. (12). This 
results from the fact that numerically 5C = 2AC for low firing rates or intracellular 
sodium concentrations (Fig. 2B). It is also interesting to note that the Attwell-Laughlin 
(2001) formula does not account for changes in the intracellular Na + concentration due 
to repetitive firing, as opposed to Eq. (12) that includes such changes through SC. In 
this sense Eq. (12) extends the phenomenological approach of these authors into the 
broader range of frequencies. 

3.2 Sodium build-up due to repetitive firing. 

When the neuron fires repeatedly, its intracellular sodium accumulates with every spike 
because Na + /K + pump is slow and cannot remove all Na + promptly (Fig. 3A,B). Since 
every spike introduces C/ n FA[Na] of sodium electric charge, and this process is fast, 
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we can approximate Eq. (6) for the Na + dynamics as 



U nF tM „ h + U„FA { N,] o p(t - U ) - jgf^, (13) 



where I Q is the sodium current associated with sodium channels and synaptic contribu- 
tion at rest, and it is given by I Q = {g Na , + q(f)fp s dg s T s V K / (4(1 - <f>)(V K - V Na ))} S(V Na - 
V ). The delta functions present in Eq. (13) represent spikes of Na + influx at times 
U. Simulation of Eq. (13) is shown in Fig. 3C, and it resembles the simulation of the 
original Eq. (6) (see Fig. 3A). 

We can further simplify Eq. (13) at the long-time limit, in which we substitute 
for the right hand side of Eq. (13) its temporal average. In particular, the temporal 
average of YLi d(t — ij) is equal to the firing rate /, and thus 



uyJM^ lo + fUnFAmo -M§!±, ( u) 

dt [Na] + 9 k 



where [Na] is the temporal average of [Na]. A comparison of the time dependence of 



[Na] from Eq. (14) with the time dependence of [Na] coming from Eq. (13) is presented 



in Fig. 3C. The equilibrium value of [Na], denoted as [Na] ot , is determined from the 



condition d[Na]/dt = 0, with the help of Eq. (12). [Na] ot) satisfies the following equation: 
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Note that the right hand side of Eq. (15) also depends on [Na] a „ through V^ a and 
SC. Thus, we can find [Na] a „ only numerically, either from Eq. (15) or from a direct 
simulation of Eq. (14). 

In Fig. 4, we compare the dependence of [Na] ot , on firing rate that comes from Eq. 
(15) with the dependence that comes from a direct numerical integration of Eq. (6). 
Overall, the formula (15) provides a relatively good approximation to the numerical 
solution, especially for low firing rates (Fig. 4). For this reason, it is used in the 
following sections to approximate metabolic expenditure and thermal changes in the 
brain tissue. 

It is interesting to note that the average [Na] a „ is weakly dependent on fiber diameter, 
especially for small firing rates (Fig. 4). The reason for this is that only the current I D 
in Eq. (15) contains synaptic contribution with d, and for low frequencies I Q is small. 

It should be kept in mind, however, that sodium fluctuations around [Na] ot , do grow 
with decrease in fiber diameter (Fig. 3B), because for thin fibers the amplitude of 
sodium influx is large (A[Na] ~ 1/d) and the relaxation time constant is short (r ~ d; 
see Appendix B). Enhanced intracellular sodium fluctuations can effectively reduce the 
concentration gradients across neuron's membrane to values close to zero, and this would 
have a devastating effect on neuron's functionality if high firing rates were maintained 
for a prolonged period of time. For example, for d = 0.11 /im corresponding to the 
smallest physically possible axons of d a = 0.06 /im (Faisal et al, 2005) the intracellular 
sodium concentration can occasionally peak to extracellular levels (145 mM; Hille 2001) 
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just for 40 Hz of repetitive firing. For d = 0.45 /iin (corresponding to d a = 0.3 /im) 
this takes place for ~ 50 Hz. These considerations suggest that very thin fibers are not 
beneficial for neuron's electrical properties. 

3.3 ATP utilization rate of Na + /K + pump, glucose metabolism, 
and firing rate in mammals. 

For pumping out Na + and pumping in K + , the Na + /K + pump uses ATP molecules. 
The average ATP utilization rate of a single neuron per its surface area S is equal to 
Ip,av/(FS). Generally, it increases with [Na] ot) and thus with the firing rate / in a 
non- linear fashion (Fig. 5A). For low frequencies / there exist a linear regime, while for 
high /, the ATP rate saturates reaching its maximum value A/F (Fig. 5A). Typical 
values of ATP utilization for frequencies in the range 1 — 10 Hz are (0.6 — 9) • 10 4 ATP 
molecules per fim 2 per second. In the linear regime, the increase of firing rate by 1 Hz 
leads to the increase of the ATP rate by about 1.7 ■ 10~ 4 /xmol/ (cm 2 sec). It is important 
to point out that in previous phenomenological models of ATP utilization rate (Attwell 
and Laughlin 2001; Lennie, 2003) only linear regime was assumed in estimations. Thus, 
the present more detailed model extends these calculations into the non-linear regime. 

The average firing rate in mammalian brains can be determined indirectly from the 
cerebral glucose utilization rates CMR 9 ;„ (expressed in mol/(cm 3 -s) ). If we assume 
that the ATP activity of the neural pumps constitutes the major contribution to the 
gray matter metabolism (Astrup et al, 1981; Erecinska and Silver, 1989; Rolfe and 
Brown, 1997; Ames 2000), then we can relate directly the ATP rate or the pump 
current to CMR giu . Since 31 ATP molecules are produced per 1 glucose molecule 
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(Rolfe and Brown, 1997), we have for the whole gray matter the following equality 
31UgCMR g i u = NI Pt a V /F, from which we obtain, using Eq. (4), that 



HMD _ 4(l-0)A[Na]j 

CMRdlu ~ 31Fd([NaL + 6>*)- (16) 



We can write CMR 9 ( U in an equivalent form, which contains neurophysiological param- 
eters explicitly. Using Eq. (15), we obtain: 



CMR glu » 4( 9 1 3Frf 0) (7 /S- + f(C + 5C){V Na - V )) . (17) 



Dependence of CMR 9itt on firing rate / is plotted in Fig. 5B, and it is practically 
the same as the dependence of I p>av on /, as the two quantities CMR 9 ; U and I p>av are 
proportional. The non-linear part of the dependence comes from the fact that 5C and 
Vjvo decrease for high frequency (via [Na] „). The CMR giu vs. / relationship (Eq. 17) 
enables us to find average frequencies for several mammalian species for which empirical 
values of CMR 9 ; U are known (Fig. 5C). In general, the average firing rates are rather 
low, from ~ 1.7 Hz for human to ~ 6.2 Hz for mouse (Table 2). Moreover, estimated 
in such a way average firing rates scale with gray matter volume with an exponent of 
—0.15 (Fig. 5C), implying that average activity in larger brains is slower than in smaller 
brains. 
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3.4 Biphasic dependence of pump power on frequency. 

The electric power generated by the Na + /K + pump in a single neuron (per surface 
area) is determined in two ways. First, from a direct numerical integration of Eq. (7) 
with time dependent voltage V and pump current I p . Second, from the approximate 
analytical formula (8) with the help of derived Eq. (15). Both methods yield similar 
results (Fig. 6A), which indicates that the approximation (8) is reliable, especially for 
low firing rates. We can rewrite Eq. (8) for the total electric power generated in the 
gray matter in a more convenient form using Eq. (4). The result is 



ATP ~ d([Na]*, + 0*) {WNa ~ VK " Vo) ■ ( 8) 



Alternatively, we can use Eq. (15) to relate the pump current to the neurophysiological 
parameters. In this way, we obtain: 



Patp 



4(1 - 4>)U g 



3d 



(3V Na - 2V K - V ) (V Na - V ) 



X 9Na,o + 



fqp s dg s T s V K 
4{l-<j))(V K -V Na ) 



+ f{C + 5C)\. 



(19) 



Note that synaptic depression via q reduces the power Patp- In Eq. (19) the first 
term in the large bracket represents a very small sodium influx at rest, the second term 
corresponds to the background dendritic synaptic activity, and the last term comes 
from Na + influx due to action potentials. The relative contribution of these 3 elements 
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to Patp depends on firing rate /. For example, for / = 1.7 Hz, corresponding to the 
estimate for human brain (Table 2), sodium influx at rest constitutes 5 %, background 
dendritic synaptic activity 15 %, and Na + influx due to action potentials yields 80 % 
of the total power. 

Fig. 6A indicates that Patp depends biphasically on firing rate /. This is a non- 
intuitive result, following from the fact that Patp is a product of two terms: I p and 
(3Vjv a — 2Vk — V ). The first of them increases monotonically with /, whereas the second 
decreases with / because Vno, and Vk depend on frequency via [Na] ot) . This biphasic 
dependence of Patp on frequency has interesting implications for thermal properties 
of brain tissue, which are discussed in the subsection (3.7). From the estimated firing 
rates for several mammalian species we can also estimate, based on Eq. (19), their 
Patp rates (Table 2). 

The sodium pump power Patp also depends inversely on the fiber diameter d (Eq. 
(19) and Fig. 6B), and proportionally on the gray matter volume U g . Thus, too thin 
fibers are metabolically expensive. 

The power generated by the sodium pumps in the gray matter can be related directly 
to the glucose cerebral metabolic rate CMR 9 ; a if we combine Eqs. (16) and (18). The 
resulting relationship is: 

Patp ~ 31FU g (3V Na - 2V K - V ) CMR glu . (20) 
Thus, the power generated scales linearly with the glucose consumption rate. However, 



23 



it should be kept in mind that the proportionality factor is not a constant, but changes 
with firing rate. 

3.5 Efficiency of the Na + /K + pump. 

Let us estimate the efficiency of the sodium pump, i.e., how much energy does it use 
for pumping out 3 Na + and pumping in 2 K + ions given an available energy from ATP 
hydrolysis. Energy from hydrolysis of 1 ATP molecule goes for performing the useful 
work of — 3e(V — Vno) + 2e(V — Vk), or equivalently 1 mole of ATP performs the work 
of F(3Vjva — 2Vr- — V ), where e is the electron charge. On the other hand, hydrolysis 
of 1 mole of ATP generates Jatp of free energy. Thus, the efficiency of the process is 
given by F(3V/v a — 2Vx — V )/Jatp- The value of Jatp depends to some extent on 
the internal chemical (ionic) state of the cell, and it has been reported to be in the 
range from 48 kJ/mol (Jansen et al, 2003) to 62 kJ/mol (Erecinska and Silver, 1989). 
This leads to the pump efficiency of 73 — 95 % (for typical resting values of voltages: 
V Na = 0.068 V, V K = -0.100 V, V = -0.067 V; Hille (2001)). Because of the high 
efficiency of the Na + /K + -ATPase, in what follows, we make an approximation in which 
we equate the heat released in the gray matter due to hydrolysis of ATP (G A tp in Eq. 
9) with the electrical power dissipated by the pump (Patp in Eq. 7). 

We can estimate the heat rate for the gray matter of human brain from Eq. (20). 
Taking U g = 680 cm 3 (Stephan et al, 1981), CMR glu = 5.7- 10~ 9 mol/(cm 3 -sec) (or 0.34 
/imol/(cm 3 -min); e.g. Clarke and Sokoloff, 1994), and for the above values of voltages 
we obtain Patp ~ 5.5 Watts. Given a possible increase of this heat value by up to 
27% due to pump efficiency, this result does not differ much from other estimates of 
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heat in gray matter. For example, Aiello and Wheeler (1995) used mass specific heat 
generation of 11.2 W/kg (based on older experimental data), which yields 7.8 Watts for 
the human gray matter. 

3.6 Scaling of cerebral blood flow with brain size. 

Cerebral blood flow CBF is important in controlling brain temperature (see the next 
subsection). The dependence of physiologically averaged CBF on brain volume can be 
found from the empirical data available in the literature. The results in Fig. 7 for 6 
mammals spanning 3 orders of magnitude in brain size show that CBF scales systemat- 
ically with brain volume as CBF = 0.018f/ 6 ~ 010 1/sec. This implies that cerebral blood 
flow decreases weakly as brains increases in size, ranging from ~ 0.020 sec -1 for mouse 
to ~ 0.009 sec -1 for human. 

3.7 Brain temperature vs. frequency and fiber diameter, and 
efficiency of brain cooling 

The spatial distribution of brain temperature T(r) is found by solving Eq. (9) with 
G atp = Patp-, and the result is given by (see Appendix C): 



T (r) = T + PATP - {aSB{T ° c ~ T ° 4) + v{Tsc ~ T °» r ^-r) 
{ } pkCmC^CBF (Kp wCw CBF)V2 



(21) 



where £ = (puCbiCBF / k) 1 ^ 2 . The parameter £ characterizes the inverse of the length 
of "transition" region from the scalp to the brain's interior where temperature is inho- 
mogeneous (Fig. 8A). For human 0.4 cm, and it is much smaller than human 
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brain radius (~ 8 cm), and thus heterogeneity is present only at the brain's edge. On 
the contrary, for mouse £ _1 ~ 0.25 cm, i.e., it is comparable with the mouse brain 
radius (~ 0.4 cm), implying that in very small brains cerebral temperature is inhomo- 
geneous in the whole volume (Fig. 8A). The scalp temperature T sc present in Eq. (21) 
is determined self-consistently (see Appendix C) from the condition T sc = T(R), and it 
decreases very weakly with brain volume (Table 2). Its value for human, T sc = 34.7°C, 
is similar to experimental values (34 — 35°C; Hensel et al, (1973)). 

Deep brain temperature T(0) is slightly larger than the blood temperature by 0.1 — 
0.2 °C for all analyzed mammals, except mouse (Table 2). This value lies in the range of 
values observed experimentally (Hayward and Baker, 1968; Nybo et al, 2002; Kiyatkin, 
2007). In general, however, T(0) is very weakly species specific. For large enough 
brains T(0) is approximately T(0) ~ T bi + Patp / (PbiCuUb r CBF) , which shows that the 
boundary temperature (and scalp cooling) becomes unimportant in this limit. From 
this formula it follows that the deep brain temperature is greater than the blood (or 
core body) temperature by the quantity proportional to the ratio of CMR S ;„ to CBF, 
in agreement with Yablonskiy et al (2000). This result suggests for example that, for a 
given animal, any local increase in the cerebral blood flow that exceeds an increase in 
glucose utilization rate leads to a lowering of brain temperature, and vice versa. 

Deep brain temperature T(0) depends significantly on the effective fiber diameter d 
(Fig. 8B,C). For very small d (corresponding to very small d a ), the temperature T(0) 
tends to diverge, which is a direct consequence of the fact that Patp also diverges for 
d — > (or equivalently d a — > 0). In this "thin fiber" limit we have a simple inverse 
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proportionality (see Eqs. (19) and (21)) between brain-body temperature difference 
and fiber thickness, i.e., T(0) — Tu ~ 1/d ~ l/(2d ). This relationship implies that 
decreasing fiber diameter twofold increases T(0)— T w by the same amount. This suggests 
again that too thin fibers are not beneficial for brain thermal equilibrium and hence its 
functioning. 

In contrast to the monotonic dependency of T(0) on d, its dependence on firing rate 
is biphasic (Fig. 8D). The origin of this non-monotonic behavior is the dependence of 
Patp on /, which is also biphasic (see Fig. 6A). This result has surprising thermody- 
namic consequences, namely too high levels of neural activity (firing rate) can lead to 
a decrease in the cerebral tissue temperature. For human brain, corresponding changes 
in temperature are rather small, at the peak about 0.5-0.6 K (Fig. 8D), which falls into 
the range of values reported experimentally (Yablonskiy et al, 2000; Kiyatkin, 2007). 

In the superficial regions brain temperature is always smaller than the blood tem- 
perature (Fig. 8A; Table 2). This fact has important consequences for the heat transfer 
in the mammalian brains. In deep brain regions cerebral blood flow plays the role of 
a coolant, whereas in the superficial regions it serves as a strong brain heater. Conse- 
quently, the net effect of blood flow is to warm up the brain tissue, i.e., Qu < (Table 
2). 

In general, the cooling mechanisms discussed in subsection (2.5) depend on brain 
temperature distribution, on scalp temperature T sc , and on brain volume Uf, r . Therefore, 
the efficiency of cerebral cooling is species specific. In particular, the relative importance 
of the two major mechanisms of heat transfer inside the brain is dependent on brain 
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size. For very small brains, heat cooling rate via conduction Q c is comparable to 
the heat warming rate via blood flow Qu (Table 2). The reason for such a strong 
warming through the cerebral blood is that it must compensate heat loss due to scalp 
convection/conduction Q cv and radiation Q r , which are much larger than the metabolic 
rate Patp- On the other hand, for large brains, the cooling rate Q c is almost twice as 
large as the warming rate Qu, because heat production due to metabolic activity Patp 
is more significant than for small brains (Table 2). On the scalp, the heat transfer rate 
is dominated by convection/conduction Q cv , as it is twice the radiation rate Q r . 

3.8 Thermal bounds on fiber diameter and length. 

Mammals are able to sustain brain temperatures up to about 42 °C without causing 
brain damage (Gordon, 1993; Kiyatkin, 2007). Above these temperatures molecular 
changes in neurons and synapses become critical and irreversible. In what follows, 
we estimate the thermal bounds on fiber diameter that allow to maintain the safe 
temperature regime. For this purpose we use Eq. (21) for r as 0, in which we put 
AT max = T(0) — T b i — 5 K, as the largest possible difference between brain and blood 
temperatures. Solving Eq. (21) for the effective fiber diameter d, we obtain 



d>d . *0--<i>)U 9 (9Na t o + f(C + 6C)) 

- mm 3 Pbl c bl U br CBFAT max _ fgp s g s T s V K U g ' 1 1 

(V N a-Vo)(3V Na -2V K -Vo) (V K -V Na ) 



where we used Eq. (19) for Patp- The minimal effective fiber diameter d min depends 
very weakly on gray matter volume U g (Fig. 9A). This suggests that brain size is not a 
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critical factor determining thermodynamic safety of the cerebral tissue. In other words, 
thermal properties of mouse and elephant brains, differing by 4 orders of magnitude in 
size, are rather similar. 

The minimal effective fiber diameter d min depends stronger and biphasically on 
firing rate / (via V^ a , Vk, and [Na] ot; ) (Fig. 9B). For very low firing rates (/ — > 0), the 
expression (22) simplifies, and we obtain in this limit: 



d _ 4(1 - <l>)g Na>0 (V Na - V )(3V Na - 2Vk - Vo)U g ^ 
3pbiCbiU br CBFAT rnax 



which yields ~ 0.3 nm (values of U g and Ut, r are for human brain). For intermediate 
values of /, for which d m i n has a maximum, the value of d m i n can be in the range 0.04 — 
0.05 /im (Fig. 9B), which corresponds to a bound on axon diameter d a = 0.02 — 0.03 jum. 
The latter value is only about 5 times larger than the membrane thickness (Koch, 1998). 
The term proportional to fq in the denominator of Eq. (22) is the background synaptic 
contribution with frequency dependent depression inside q. Due to this depression the 
synaptic contribution is bounded from above and consequently does not diverge as a 
function of /. Thus, synaptic depression not only reduces the power Patp (see. Eq. 
19), but it also makes d min finite. Without depression, the bound d min would diverge 
for / ~ 450 Hz. 

Average value of the empirical fiber diameter d in the gray matter is brain size inde- 
pendent and is about 0.45 /im (harmonic mean of average axon and dendrite diameters 
in mouse, as defined in Sec. 2.2, Eq. (3)). This value is about 10 times larger than 
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the largest value of d m i n . Because of this, mammalian brains operating under normal 
physiological conditions are rather safe from excessive overheating that would cause 
brain damage, for all ranges of frequency. The situation could be more tricky in hot en- 
vironments in which AT max were severely reduced, due to increase in body temperature 
and hence cerebral blood temperature. The margins of thermal safety, could be also 
compromised in pathological conditions. For example, abnormalities associated with 
strongly reduced cerebral blood flow CBF or compromised synaptic depression, could 
significantly increase the bound d m i n up to the lower neuroanatomical values of d. 

The thermal lower bound on fiber diameter also determines the upper bound on 
fiber length per neuron I, as is evident from Eq. (4). Taking the neuron's surface 
area S = nld, we get / < 4(1 — (f>)U g /(7iNd^ nin ), where d min = 0.01 — 0.05 jum. For 
human brain with U g = 680 cm 3 (Stephan et al, 1981) and N = 2 • 10 10 (Haug, 1987; 
Braendgaard et al, 1990) this yields I < 12—300 meters. For mouse brain with U g = 0.11 
cm 3 and N = 1.6 ■ 10 7 (Braitenberg and Schuz, 1998) we get I < 2.4 — 60 meters. 
Neuroanatomical data for mouse gray matter indicate that the average fiber length per 
neuron is 1.5 — 4.5 cm (Braitenberg and Schuz, 1998), which is 2-3 orders of magnitude 
below its upper thermal bound. A similar conclusion holds for the human brain. 

Taken together all these results suggest that thermodynamics does not restrict neu- 
roanatomical parameters in any dramatic way, because they are far away from their 
thermal bounds. 
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4 Discussion 



This article investigates thermodynamic properties of brain tissue and corresponding 
physical limits on neural anatomy caused by heat balance in the brain. It is found 
that, in general, the lower and upper limits on fiber dimensions are unattainable for 
normal values of physiological parameters such as cerebral blood flow and maximal 
Na + /K + -ATPase. This suggests that real mammalian brains either keep these physio- 
logical parameters in the "proper" range or scale appropriately average fiber diameter 
and length to maintain wide margins of thermodynamic safety. Such wide margins 
presumably enable the brain to avoid overheating during enhanced cerebral activities, 
abnormal states (e.g. epileptic seizures), and in hot environments. 

The conclusions of this paper are based on calculating neural metabolic power and 
spatial distribution of heat dissipated. Since majority of metabolic energy in neurons 
goes to pumping out sodium ions (Ames, 2000; Astrup et al, 1981; Erecinska and Silver, 
1989; Rolfe and Brown, 1997), the first step in determining neural power was to find 
Na + influx during an action potential. An explicit analytical formula for this Na + 
influx was derived (Eq. 12) and was validated numerically. This analytical formula 
yields similar results for an isolated action potential as the one of phenomenological 
character by Attwell and Laughlin (2001) (20 — 28% discrepancy between the two can 
be attributed to the different assumptions regarding the magnitude of depolarizations 
in axons and dendrites). However, for repetitive firing the derived here Eq. (12) extends 
the Attwell and Laughlin (2001) formula, because it accounts for intracellular sodium 
accumulation and corresponding changes in membrane electrical properties expressed 
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by SC (correction to Na + influx due to prolonged Na + channels activation). Moreover, 
the sodium build-up is also important for the accurate determination of ATP/glucose 
utilization rates for high frequencies. In this respect, Eqs. (16) and (17) extend the 
Attwell and Laughlin (2001), and Lennie (2003) linear calculations to the non-linear 
regime. An additional novelty in the present approach is inclusion of synaptic depression 
in the transmission probability. 

The issue of brain cooling is not a classic problem of volume to surface ratio, as it is 
the case with natural, non-designed, physical objects. Instead, the brain cooling could 
be compared to the cooling of combustion heat engine, which receives a liquid coolant. 
In the brain the role of the coolant is played by the cerebral blood, but only in the deep 
region because there blood has a slightly lower temperature than the brain tissue. In 
the superficial regions brain tissue has a smaller temperature than the cerebral blood, 
and there blood warms up the brain. The fact that the deep brain temperature depends 
weakly on brain volume (Table 2), implies that brain size is not a major determinant of 
thermal responses. This in turn implies that the thermodynamics of heat balance does 
not restrict the brain size in any significant way, suggesting that, in principle, brains 
could be heavier than 5 kg (the largest known brain). 

The interesting result is that the power generated by the sodium pump depends 
biphasically on neural frequency of firing, and inversely on the effective fiber diameter 
(Fig. 6). As a consequence of this, brain temperature T(0) can depend biphasically 
on frequency as well, if cerebral blood flow does not change with frequency (Fig. 8D). 
Thus T(0) increases with frequency but only up to a certain point above which it 
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slightly decreases with further increase in frequency. The increase of brain temperature 
in response to activation (higher firing rate) is an expected behavior (Kiyatkin, 2007). 
However, the decrease of T(0) for very large stimulation (very high firing rates) is an 
unexpected effect that can explain some experimental results in which stimulation of 
certain brain regions led to lowering of local temperature (McElligott and Melzack, 
1967; Yablonskiy et al, 2000). The standard explanation for the decrease of T(0) upon 
stimulation requires that cerebral blood flow must increase more than Patp increases 
(Yablonskiy et al, 2000). The biphasic relationship between temperature and frequency, 
obtained in this paper (Fig. 8D), offers an alternative explanation, namely, that brain 
temperature can decrease upon vigorous (or prolonged) stimulation even when cerebral 
blood flow does not change. 

The negative correlation between the power generated in the brain and the effective 
wire diameter (Fig. 6B) suggests that thin fibers are energy expensive. The effective 
fiber diameter d is defined as a harmonic mean of axon d a and dendrite da diameters (Sec. 
2.2), with the assumption that volumes of axons and dendrites are equal (Braitenberg 
and Sciiz, 1998; Chklovskii et al, 2002). In general, d and d a are positively correlated 
and c? — ^ if <i a — > 0. Thin fibers also promote large fluctuations of intracellular 
sodium concentration in response to stochastic action potentials (Fig. 3B), which can 
be disadvantageous for electrical properties of neurons. 

Another negative consequence of having too thin fibers is that cerebral tissue tem- 
perature increases inversely with fiber diameter (Fig. 8C), i.e., T(0) — T b i ~ 1/d. Thus, 
brain regions reach in very thin wire can heat up excessively. As an example, the tem- 



33 



perature of the thinnest known axons with diameter 0.1 jum (Faisal et al, 2005, i.e. 
d — 0.18 /im), relative to blood temperature, should be about 2.5 times larger than the 
corresponding relative temperature of 0.3 /im axons (mouse cerebral cortex; Braiten- 
berg and Schuz, 1998; corresponding to d — 0.45 /im), assuming both axons have the 
same firing rate. Thus, if typical relative temperatures above blood temperature are 
0.2 — 0.3 °C for 0.3 mum axons (Hayward and Baker, 1968), then for 0.1 /im axons 
the corresponding relative temperatures would be 0.5 — 0.8 °C. This result suggests an 
explanation why in the peripheral nervous system sensory fibers responsible for high- 
threshold heat sensation (so-called C-fibers) are much thinner (0.1 — 1.5 /im; Kandel 
et al, 1991) than other sensory fibers (A- type or B-type with diameters 1 — 5 /im). It 
might be that they warm up easier, although these fibers also serve other functions 
(Craig, 2003). 

This study finds that the lower thermal bound on the effective fiber diameter d m i n 
is strongly frequency dependent, but it is finite due to synaptic depression. Values of 
dmin are in the range 0.3 — 50 nm (or 0.0003 — 0.05 /im; see Fig. 9B). These values can 
be translated to the corresponding thermal bounds on axon diameter d a , which are in 
the range 0.0002 — 0.026 /im (assuming that the dendrite diameter d^ = 0.9 /im and it 
does not change). The average value of cortical axon diameter is d a = 0.3 /im, and thus 
it is 12-1500 times larger than these limits, therefore, on average, mammalian brains 
operate in the safe thermal zone. However, it should be also kept in mind that axon 
diameter displays some variability, and the thinnest axons can reach 0.1 /im (Faisal 
et al, 2005). This value is still about 4 times larger than the maximal value of the 
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lower thermal limit on d a (which is 0.026 /im), suggesting that thermal limits are not 
an immediate constraint on fiber size. It is also interesting to note that the estimated 
maximal lower thermal bound on the diameter of brain axonal wiring (0.026 /im) is 
smaller than corresponding bounds imposed by structural constraints and noise, which 
are respectively 0.06 /zm and 0.1 /im (Faisal et al, 2005). 

The theoretical bound on axon diameter implies a corresponding lower bound on 
the speed of signal propagation in the gray matter. Let us estimate the upper limit on 
temporal delays. For unmyelinated fibers (prevalent in gray matter) velocity of signal 
propagation is proportional to the square root of the fiber diameter (Hodgkin, 1954). 
Experimentally, for axons with d — 1 /jm the propagation velocity is 2.3 mm/msec 
(Koch, 1998), which implies that for axons with the boundary thickness of 0.026 /im, 
we have velocities 0.37 mm/msec. This gives for the maximal known extent of axons 
in the gray matter of 9-10 mm (for macaque monkey visual cortex; Amir et al, (1993)) 
the upper limit on delays in the range 24-27 msec. Thus, apparently, thermodynamics 
of heat balance in the gray matter does not tolerate temporal axonal delays longer than 
~ 0.03 sec, which is a stringent constraint. 

The problem of finding bounds on fiber dimensions is similar in spirit to the ap- 
proaches of "wire minimization" in the brain (Cherniak 1995; Murre and Sturdy 1995; 
Karbowski 2001, 2003; Chklovskii et al, 2002). It is hypothesized that this principle 
governs the organization of the mammalian nervous system at different scales (Murre 
and Sturdy 1995; Prothero, 1997; Kaas 2000; Karbowski 2001, 2003), because it of- 
fers energy savings associated with ionic membrane transport, as well as reduction in 
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temporal delays in neural communication. The heat balance limits on fiber diameter, 
found in this paper, suggest corresponding upper thermal bounds on the density of wire 
packing in the brain in the range 20 — 100 fibers per /iin. These thermal bounds on wire 
density are a factor of 2-4 higher than the corresponding upper bounds coming from 
structural and noise considerations (Faisal, 2005). 

The estimated firing rates in mammals are in the range from 1.7 Hz for human to 6.2 
Hz for mouse, and they scale systematically with brain size, with the exponent —0.15 
(Fig. 5C). The estimate for rat (5 Hz) is very close to that assumed by Attwell and 
Laughlin (2001), i.e. 4 Hz, which was based on weighted average of values observed 
experimentally. Also the estimate for cat (4.5 Hz) is reasonably close to that reported 
experimentally for visual cortex (2.5-4.0 Hz; see Baddeley et al, 1997). However, in vivo 
values of firing rates are stimulus dependent and in many species are largely unknown. 
The scaling result of firing rates is qualitatively consistent with the experimental data on 
avian brains, which show an allometric decay of firing rates in the peripheral nervous 
system with brain/body mass (Hempleman et al, 2005). It is also consistent with a 
prediction coming from a recent analysis of empirical data on brain metabolic scaling in 
mammals (Karbowski, 2007), where it was suggested that an average firing rate should 
decrease for bigger brains to account for a negative allometric exponent of specific 
metabolic rate, which was also —0.15. This conclusion is analogous to the general trend 
in mammals, in which physiological processes tend to slow down with an increase in 
body size (Schmidt-Nielsen 1984). 

The estimate of the heat generated in the gray matter may have some margins of 
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error. First, the assumption of the equipotential neuron is only an approximation, 
because it does not include explicitly the spatio-temporal effects associated with action 
potential propagation and back-propagation. However, it is estimated in Appendix A 
that the correction from this effect is rather small of the order of 16%, i.e., the actual 
Na + influx can be larger by this amount. Second, because the efficiency of the pump 
is a little less than 100% (73 — 95 %; Sec. 3.5), the actual heat dissipated in the gray 
matter may be larger than Patp by 5 — 27 %. Therefore the total error from these two 
effects on the cerebral heat can theoretically reach 43%. We should also remember that 
some heat coming from the glucose to ATP conversion was neglected (i.e., it is assumed 
that glycolysis is 100% efficient). However, the comparison of the heat estimates for the 
human gray matter indicates that in fact the total heat error cannot be larger that 30% 
(Sec. 3.5). There are also other sources of error that affect the temperature distribution 
in the brain, such as the neglect of cooling by scalp perspiration, and variability of 
environmental T Q temperature, which can affect convection/conduction and radiation 
cooling rates. However, these contributions should not have any dramatic influence on 
the thermal properties of brain tissue. Moreover, neither of these sources of error or 
their combination seem to affect the main conclusion of this paper, namely, that neural 
anatomy (fiber diameter and its length) does not approach its thermodynamic limits. 
To reach these limits the ratio of the maximal sodium pump current to the cerebral 
blood flow would have to be at least 10-100 times larger. 

The formulas (Eqs. 17-21) in this paper for the CMR g ; u , power dissipated in the 
gray matter, and brain temperature, may have practical use. These formulas as well as 
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their possible future extensions can be used for assessing neural activity (firing rates) 
based on changes in temperature, CMR g i u , and CBF. Moreover, because this study 
offers a direct relationship between neural activity and neuroanatomical parameters, 
it may be useful in quantitative studies of the interplay between brain development, 
evolution, and metabolism (Purves, 1988; Striedter, 2005). 
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Appendix A: Sodium influx during an action potential. 

The duration of a typical action potential can be divided into two phases (Fig. 
2A). During the first phase Na + conductance g^ a rises almost instantenously to its 
maximal value g Na and voltage V increases to its peak value V^ a . The second phase is 
characterized by decline in values of gNa (to zero) and V (to values < 0). The total Na + 
influx during an action potential is A[Na] G = A[Na]^ + A[Na]^ 2 \ where superscripts 
(1) and (2) refer to the first and second phase, respectively. 

During the first phase, the intracellular potassium concentration practically does 
not change, because K + channels are activated with a delay. Thus, by straightforward 
integration of Eq. (2), we obtain the Na + influx during this phase as: 
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A[Na]« 



4C 
Jd 



(V Na -V ), 



(24) 



where the contribution proportional to ql was neglected, since it is much smaller. 

Sodium influx during the second phase can be computed from the sodium dynamics 
of Eq. (6). Most of the time during this phase, the term proportional to gNa is much 
larger than the remaining two terms in Eq. (6). Thus, we can write 



where r is the duration of the second phase and r c 0.4 msec (from simulations). To 
simplify calculations, it is assumed that Na + conductance and voltage depend on time in 
the following way: g Na (t) ^g Na (l-t/r ) and V(t) « V Na + (0.QV K -V Na )(t/T o ) z . These 
forms assure that for t = 0, we have gNa(fy ~ ~9~Na an d V(0) ~ Vno,, and for t = t d we 
have gNa(T ) and V(t ) ks O.QVk- The latter value comes from simulations. After 
performing integral in Eq. (25) we obtain: 




(25) 



A[Na]i 2 ) 



4^ a r o (W Q -0.6V^) 
Fd(z + l)(z + 2) 



(26) 



The fits in Fig. 4 are made for z = 2.5. The total Na + influx during an action potential 
is A[Na]W + A[Na](, 2) , and is given by Eq. (12) in the main text. 
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It is also interesting to check the magnitude of correction coming from the fact that 
real neurons are not equipotential and action potentials propagate and back-propagate 
with a finite velocity. Let ( be the spatial constant characterizing membrane poten- 
tial homogeneity, and c the velocity of action potential propagation (along axon) and 
back-propagation (along dendrites). A simple way to account for the spatio-temporal 
dependence of Na + conductance and voltage in Eq. (25) is to make the following rescal- 
ings: 

g Na (x,t) ^g Na (l - (t - x/c)/t ) exp(-(c£ - x)/()H(ct - x)H(x + cr Q - ct) 
and 

V(x,t) » {V Na (l -(t- x/cY/t*) exp(-(ct - x)/0 + 0.6V^(t - x / c) z / r z ] H (ct - 
x)H(x + cr — ct), 

where the factor exp(— (ct — x)/() denotes a traveling wave of excitation along axon 
or dendrite with the spatial spread (, and the function H(y) is the standard Heaviside 
function equal to 1 for y > and equal to for y < 0. These rescalings assure that 
points along axon and dendrite separated by x from the soma (action potential initiation 
zone) receive an excitation after the delay time x/c. With these modifications, after 
some algebra, the corrected result for the A[Na]^ is 

[ Jo Fd(z+l)(z + 2) \ +€ {)' { ] 

where e is given by e = {v n ^-q\v k ) 

{(z + 2)V Na /6 - (2V Na - 0.6V K )/(z + 3)}. The 
term proportional to cr /( represents the first order correction and it is clear that the 

40 



equipotential approximation corresponds to the case when ( = oo or ct q /( <C 1. For 
physiological values of parameters: Vno, = 60 mV, Vk = —100 mV, c = 2.3 mm/msec, 
( — 1 mm (Koch, 1998), with r D = 0.4 msec and z = 2.5, we obtain ecr /( w 0.22. This 
implies that the equipotential assumption slightly underestimates the sodium influx 
during the second phase of sodium activation. Overall, this correction to the total Na + 
influx A[Na] D is small, ~ 16%, for low firing rates for which 5C ps 2.4C, and it gets even 
smaller for high firing rates. A detailed numerical treatment of the influence of action 
potential velocity on metabolic rate of the squid giant axon is presented in (Crotty et 
al, 2006). 

Appendix B: Relaxation time constant for Na + dynamics. 

In this Appendix the relaxation process of voltage and sodium to their equilibrium 
values is analyzed following a single action potential for a neuron that prior to that has 
been at rest for a long time. Computations below are performed in a very late phase of an 
action potential when the conductances g^ a , 9k, and the voltage V returned essentially 
to their resting values (or are close to them, i.e. dV/dt ps 0). During this phase we can 
use a linear approximation on V and [Na] in Eq. (6). That is, we can expand V and [Na] 
around their resting values V and [Na] G as: V ~ V Q + AV, [Na] ~ [Na] + A[Na], where 
A[Na]/[Na] 1, and neglect higher order terms. In this approximation the pump 
current I p can be written as I p ps l V)Q + AA[Na], where I Pt0 = AS[Na]J/([Na]J + 8 k ), 
and A = kA9 k [Na\^~ 1 / \6 k + [Na]^) 2 . The corresponding changes in VVa and Vk are 
AV Na » -(i?T/F)A[Na]/[Na] and AV K » - (RT/F) A [K]/[K] = (i?T/F)A[Na]/[K] . 
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Next, using the facts that gNa/gi «C 1 and gx/gi <C 1, we can solve Eq. (6) analytically. 
The equations governing A[Na] and AV relaxations are given by: 



A[Na] = - 



3A + 



F[Na] 



RT 




AV(t) 



{A + (RT/F)(g Na /[N a ] - g K /[K})} AjNaj 
{gL + qfp s g s T s d/ (4(1-0))} 



(28) 



The term on the right hand side of A[Na] dynamics that is proportional to RT/ (F[Na] G ) 
is much smaller than 3A, and thus it can be neglected. This implies that changes in 
Vnu and Vk due to sodium influx for an isolated action potential are very small. This 
leads to a simple exponential decay of A[Na] as A[Na] xs A[Na] G e~*/ r with the time 
constant r xs Fd/(12A), where A[Na] G is Na + influx during the action potential given 
by Eq. (12). Thus, relaxation time constant is short for thin fibers, and it increases 
proportionally with fiber diameter. As an example, for d = 0.45 /im we obtain r fa 5 
sec (for k = 3 and [Na] G = 12.0), which is of the right order of magnitude (Abercrombie 
and Weer, 1978; Nakao and Gadsby, 1989). For d — 0.18 /im the time constant r is 
about 2 sec. 

Appendix C: Solution of the thermal balance equation. 

The steady-state limit of Eq. (9) can be rewritten in the form: 



d 2 T p«c«CBF 



(T-T ), 



(29) 



dr 2 k 
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where T Q = T hl + GATp/(PbiCbiCBFUb r ). In computations, we assume Gatp ~ Patp 
because of the high efficiency of Na + /K + -ATPase (Sec. 3.5). We look for the solution 
of Eq. (29) in the form: T(r) = T + aexp(£r). Substituting this form in Eq. (29), 
we obtain £ = ±(p M c^CBF//t) 1 / 2 , and we take only the positive solution as the one 
corresponding to the physical situation. From the boundary condition in Eq. (10) we 
obtain the coefficient a = - {a SB (T^ c - T G 4 ) + r](T sc - T Q )} e~^/«). Combining these 
results we find Eq. (21) in the main text. The scalp temperature T sc is determined 
numerically from the condition T sc = T(R), for all considered species. 
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Figure Captions 
Fig. 1 

Heat transfer in the brain. Brain is represented as half of the ball with three concentric 
layers representing (from the brain outside): cerebrospinal fluid, skull, and scalp. The 
heat generated in the brain is removed by cerebral blood flow (small circles) , conduction 
through the brain and three layers (dashed arrows), and scalp convection/conduction 
and radiation (dashed circular line). 

Fig. 2 

Sodium influx during an action potential. (A) Temporal dependence of Na + (solid line) 
and K + (dashed line) conductances and voltage (dashed-dotted line) during the action 
potential. (B) The correction 5C as a function of intracellular sodium concentration. 
(C) Dependence of sodium influx A[Na] Q on the effective fiber diameter d coming from a 
numerical integration of Eq. (6). The least square fit (solid line) yields the relationship: 
A[Na] D = 2.72- 10~ 5 d _a974 mM, where d is in cm. This numerical fit practically confirms 
the theoretical dependence A[Na] D ~ d^ 1 in Eq. (12). 

Fig. 3 

Dependence of the intracellular sodium concentration on time for a neuron firing re- 
peatedly. (A) Results for d = 0.45 /im, i.e., the harmonic mean of the mouse axon 
diameter d a = 0.3 fxm and dendrite diameter d d = 0.9 /im (see Sec. 2.2). (B) Results 
for d = 0.11 /xm, i.e., for the smallest molecularly possible axons with d a = 0.06 fim. 
Note that the range of variability in (B) is much larger, which reflexs higher amplitudes 
of Na + influx and its faster relaxation for thin fibers. Panels (A) and (B) come from 
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simulations of Eq. (6). (C) Simulation results for d = 0.45 jum of Eq. (13) (solid 
lines) and its approximation Eq. (14) (dashed and dashed-dotted lines). For all plots 
in (A)-(C) k = 3. 

Fig. 4 

Intracellular sodium concentration as a function of firing rate. The point (diamonds) 
coming from a direct numerical integration of Eq. (6) are approximated well by the 
theoretical formula (15) (solid line). The maximal discrepancy for high frequency is 
about 15 — 20 %. The average [Na] a „ does not depend much on fiber diameter (panels 
A and C) nor on the Hill coefficient k (panels A and B). 

Fig. 5 

ATP and glucose utilization rates, and scaling of firing rate with brain size. (A) ATP 
rate increases non-linearly with frequency / and saturates for high values of /. (B) 
CMR g i u also increases non-linearly with firing rate, similar to the dependence of ATP on 
/. (C) Scaling of the estimated firing rates in mammals with brain size. The least square 
fit to the data points yields f est = 4.79C/-°- 15 Hz (R 2 = 0.89, p = 0.0015). Estimates 
of frequency are based on the empirical data of CMR ff ; u , which are (in /imol/cm 3 min): 
1.07 for mouse, 0.90 for rat, 0.83 for rabbit, 0.81 for cat, 0.47 for rhesus monkey, 0.46 
for baboon, 0.34 for human. These data were taken from the data gathered in the 
supplementary information of Karbowski (2007). For all panels: d = 0.45 /im, k — 3. 

Fig. 6 

Dependence of the sodium pump power Patp on frequency and the effective fiber di- 
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ameter. (A) Biphasic dependence of Patp (per neuron per surface area) on frequency. 
Analytical formula (19) (solid line) leads to similar results as a direct numerical integra- 
tion of Eq. (6) (diamonds). Results are for d = 0.45 /im. (B) Patp depends inversely 
on fiber diameter d. For extremely thin fibers Patp tends to diverge. Results are for 
U g = 680 cm 3 and / = 1 Hz, corresponding to human brain. In both panels A and B, 
the Hill coefficient k — 3. 

Fig. 7 

Log-Log plot of the empirical dependence of the cerebral blood flow rate CBF on brain 
volume Ubr for several mammals. The allometric dependence has the following form: 
CBF = 0.018f/ 6 ; 010 1/sec, where U br is expressed in cm 3 (B 2 = 0.976, p = 0.0002). The 
CBF empirical data are as follows: 20.2 • 10~ 3 sec -1 for mouse (Frietsch et al, 2007); 
17.7 • 10" 3 sec" 1 for rat (Waschke et al, 1993; Linde et al, 1999); 14.3 • 10~ 3 sec" 1 for 
rabbit (Busija, 1984); 12.7- 10~ 3 sec' 1 for dog (Marcus and Heistad, 1979); 10.7- 10~ 3 
sec -1 for sheep (Koehler et al, 1985); 9.2 • 10~ 3 sec" 1 for human (Madsen et al, 1991). 
Brain volumes were taken from Stephan et al (1981) and Karbowski (2007). 

Fig. 8 

Brain temperature as a function of brain size, fiber diameter, and firing rate. (A) The 
spatial distribution of brain temperature T(r). For very small brains T(r) is smaller 
than the arterior blood temperature T u (dotted line). (B), (C) T(0) strongly increases 
for very thin fibers (U g = 680 cm 3 corresponding to human brain for panel B and C; 
and / = 1.7 Hz for panel B). (D) Biphasic dependence of T(0) on firing rate (d = 0.45 
/im, k — 3). Solid line corresponds to human and dashed line to mouse. 
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Fig. 9 

Lower bound on the effective fiber diameter d as a function of brain size and firing rate. 
(A) d min depends very weakly on gray matter volume. (B) The lower bound on d min 
depends biphasically on firing rate (U g = 680 cm 3 ). Weak synaptic depression (7 = 0.1; 
dashed line) leads to a similar dependence as moderate synaptic depression (7 = 0.5; 
solid line). The lower thermal bound on axon diameter d a is approximately two times 
smaller than d m i n , and therefore it is maximally 0.026 jum. 
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Table 1: Parameters used in the article. 



Parameter 


Value 


Units 


Reference/ Comment 


n at 
yNa,o 


2.9 ■ 10~ 7 


(Qcm 2 )^ 1 

\ u i. V XXX 1 


p\t, — fiT xy\\1 TVft nil- A/Til ps innHpl 1 1 QQ1 I 




1.2 • 10~ 6 


(^Icm 2 )' 1 

\ u Lv ill 1 


— fiT yt\\1 Tr^iih-TVTilpt; innHpl i1QQ1 i 


iJNa 


fl 1 

\J. J. 


(Orm 2 )^ 1 


Tranh-Milp^ moHpl (1 QQ1 ^ 

licLUIJ IVlllco IxlUU-Cl I XC/c/J.) 


y k 


0.085 


(ilcm 2 )' 1 

\ u u V XXX I 


Traub-Miles model (1991) 

xi dLiu ivxiivju xxx v.y v.x v x 1 Xi^i/x j 


n r 


10~ 4 


(rtcm 2 )- 1 


Tranb-Milps mnHpl f 19911 


Or. 

its 


0.3 ■ 10~ 9 




Yoshimura et al (1999) 


fin 


0.17 


unitless 


Vokushev et al (2004) 

V W X ■LXl-IXXV.j V V_j U CA X 1 — \J V 7 XI 1 








Markram et al ('1998") 

±\±C\il. xVX OjIII C u CXI I X JJU ( 








Yoshiinnra pt h1 f 1 QQQ ) 

X V/iJXllXlXl L X CXi V_. U CXiX l J- kJ kJ kJ 1 


v 

v s 


0.000 


v 


Koch (1998) 


c 


10 -6 


F/cm 2 


ffenerallv accented 

CL VjllVl CXiXX y CXiV_,\_, V_j 1J U Vj v_x 


1 s 


2.2 • 10~ 3 


sec 


Koch (1998) 


T 7 


0.5 


sec 


Davan and Abbott (2001) 


'Y 
l 


0.5 


unitless 

LXXXX UXV-JkJkJ 


Davan and Abbott f2001) 

J s CXi y CXiXX CXiXX V_X X i V-J \-J v_y U U V ^ \J \J x^ / 


A 


2 • 10~ 6 


A/cm 2 


Abercrombie and De Weer fl978) 








Nakao and Gadsby (1989) 


r S 


5 • 10 11 


cm -3 


Braitenberg and Schuz (1998) 








Craig (1967) 








DeFelipe et al (2002) 




1/3 


unitless 


Braitenberg and Schuz (1998) 


F 


96485.3 


C/mol 


Faraday constant 


K 


5 • 10- 3 


W/(cm K) 


van Leeuwen et al (2000) 


&SB 


5.67- 10~ 12 


W/(cm 2 K 4 ) 


Stowe (1984) 


V 


1.2 • 10~ 3 


W/(cm K) 


Nelson and Nunneley (1998) 








van Leeuwen et al (2000) 


Cbl 


3.8- 10 3 


J/(kg K) 


generally accepted 


Pbl 


1.06 • 10- 3 


kg/ cm 3 


generally accepted 


T u 


309.8 


K 


cerebral blood (body) temperature 


T 


293.2 


K 


room temperature (20 °C) 



55 



Table 2: Estimated neural activities and thermal brain properties in several 
mammals. 



Species 


/(Hz) 


Patp (W) 


T(0) (°C) 


Tsc (°C) 


Qu (W) 


Q c (W) 


Qcv (W) 


Qr (W) 


Mouse 


6.18 


0.003 


36.57 


35.5 


-0.024 


0.027 


0.017 


0.009 


Rat 


5.03 


0.008 


36.69 


35.4 


-0.058 


0.066 


0.043 


0.022 


Rabbit 


4.59 


0.054 


36.81 


35.3 


-0.199 


0.25 


0.167 


0.086 


Cat 


4.47 


0.27 


36.85 


35.2 


-0.51 


0.78 


0.51 


0.26 


Macaque 


2.38 


0.53 


36.76 


35.0 


-1.25 


1.78 


1.18 


0.60 


Baboon 


2.33 


0.84 


36.76 


35.0 


-1.65 


2.49 


1.64 


0.84 


Human 


1.68 


5.41 


36.73 


34.7 


-6.02 


11.43 


7.55 


3.88 



In estimations it was assumed that blood temperature Tu = 36.6 °C, and 
it is brain size independent. The following values of gray matter volume U g 
were used: 0.11 cm 3 for mouse; 0.42 cm 3 for rat; 3.0 cm 3 for rabbit; 15.2 cm 3 
for cat; 50.0 cm 3 for macaque monkey; 80.0 cm 3 for baboon; and 680.0 cm 3 
for human. 
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Figure 8B 
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Figure 8C 
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